What Does It Take to Stabilize Gravitational Clustering? 
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ABSTRACT 

An analytical understanding of the strongly nonlinear regime of gravitational col- 
lapse has been difficult to achieve. The only insight has been the stable clustering hy- 
pothesis, which assumes that the number of neighbors for objects averaged over small 
length scales is constant in time. Our recently proposed analytic halo model for N- 
point correlation functions now provides a tool for calculating gravitational clustering 
properties in the strongly nonlinear regime. This model also provides a new physical 
framework for an independent evaluation of the validity of the stable clustering hypoth- 
esis. We derive the asymptotic nonlinear behavior of the iV-point correlation functions 
and pairwise peculiar velocities in terms of dark matter halo properties. We show that 
these statistics exhibit stable clustering only when the halo mass function and halo 
density profile obey specific relations. The long-cherished stable clustering hypothesis 
therefore is not necessarily realized in practice. 

Subject headings: cosmology : theory - dark matter - large-scale structure of universe 

1. Introduction 

Gravitational clustering is the fundamental physical process responsible for the formation 
and evolution of structure in the universe. Galaxy-hosting dark matter halos are the products 
of the strongly nonlinear stage of this process, but a detailed understanding of this important 
regime has mostly eluded us because few numerical simulations have the dynamic range to explore 
comfortably such a small length scale (e.g., Moore et al. 1999; Bullock et al. 1999). The one analytic 
handle has been the stable clustering assumption (Peebles 1974; Davis & Peebles 1977; Peebles 
1980), which has allowed us to extrapolate into nonlinear regimes beyond the reach of numerical 
simulations. Predictions of the stable clustering hypothesis are: (1) the two-point correlation 
function of the density field £(r) and the power spectrum P(k) are power laws, with £(r) cx r -7 
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and A(k) = A7rk 3 P(k) oc k J , where 7 = (9 + 3n)/(5 + n) for a primordial spectral index n; (2) the 
higher order iV-point correlation functions ^iv(r) scale as £jv(r) oc r _7JV , with 7^ = (iV — 1)7, or 
6v oc f^" 1 ; and (3) the pairwise peculiar velocity exactly cancels the Hubble expansion on small 
scales, v/Hr = —1, so that bound, high-density halos maintain a fixed physical size. Some of these 
predictions have agreed well enough with numerical simulation results (Jain 1997 and references 
therein) to have obtained general acceptance. 

In two recent papers we presented evidence from high-resolution iV-body simulations that 
gravitational clustering does not necessarily follow the scaling £3 oc £ 2 required by the stability 
condition for the two- and three-point functions £ and £3 (Ma & Fry 2000a). For a deeper un- 
derstanding of the strongly nonlinear regime beyond simulations, we proceeded to construct an 
analytic halo model (Ma & Fry 2000b) in which mass is distributed in spherical dark matter ha- 
los with phenomenological mass distribution functions, density profiles, and halo-halo correlations. 
We showed that this halo model can reproduce analytically the two- and three-point correlation 
functions measured in numerical simulations, and that it also makes useful predictions beyond the 
limited range of validity of simulations. 

In this Letter we investigate the extent to which the analytic halo model is consistent with 
stable gravitational clustering in the strongly nonlinear regime. In particular, we derive the halo 
model predictions for the asymptotic nonlinear behavior of the A-point correlation functions (§2) 
and the pairwise peculiar velocities (§3) in terms of dark matter halo properties. We then obtain 
the conditions that must be satisfied by the halo mass function and density profile in order to 
reproduce results of the stable clustering hypothesis. 



The analytic halo model of Ma & Fry (2000b) for the iV-point correlation functions of the 
mass density takes three inputs: the halo mass distribution function dn/dM; the halo density profile 
p{r)/p = 5 u(r/R s ) where R s is a scale radius and 5 is the normalization; and halo-halo correlations. 
In this model, the two-point function contains contributions from particle pairs residing in a single 
halo and in two separate halos, £(r) = £i/i(r) + ^2h{ r ) , where the subscripts "l/i" and "2/i" denote 
the respective "1-halo" and "2-halo" contributions (see also Seljak 2000). Similarly, we write the 
bispectrum B(k\,k2,ks) as three separate terms, representing contributions from particle triplets 
that reside in a single halo and in two and three distinct halos. As shown in Ma & Fry (2000b), the 
dominant contribution in the nonlinear regime on small length scales is from the "1-halo" terms for 
particles that reside in the same halos. This makes intuitive sense because closely spaced particles 
are most likely to be found in the same halo. These "1-halo" terms for the two- and Appoint 
statistics are 
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H N Mn,---,r N ) = J dM^Rl~5 N X N ( n /R s ,---,r N /R s ), (1) 
P NAh {k 1 ,...,k N ) = J dM-^iRtJuihRs)} ••• [R 3 Ju(k N R s )}. 

where \n(xi, ■ ■ ■ , %n) = / d 3 y u(y) u(\y + x± — X2D ■ ■ ■ u(\y + X\ — x^\)\ u\s the Fourier transform 
of u(x); the iV-point Pjy is defined for ^ ki = 0, and P3 = -B is the bispectrum. 

For our goal of studying the iV-point functions in the strongly nonlinear regime, we need to 
consider only the halo mass function and halo density profile and not halo-halo correlations because 
the latter make negligible contributions on small length scales and are therefore irrelevant for stable 
clustering (see Ma k, Fry 2000b for the full expressions for P and B). For the halo mass function, 
we write it as 

-^ocM-^e- 2 / 2 , u = -^fr, (2) 
dM a(M) v ' 

where a parameterizes the uncertainty in the logarithmic slope at the low-mass end, S c ss 1.68 
characterizes the linear overdensity at the onset of gravitational collapse, and a(M) is the linear 
rms mass fluctuations in spheres of radius R, where M = AirpR 3, /3. Values in use for a vary from 
1 (Press & Schechter 1974) to about 0.4 (Sheth & Tormen 1999; Jenkins et al. 2000). For the halo 
density profile, we assume a "universal" shape, p/p = 5u(r/R s ), where u = l/[x(l+a:) 2 ] is proposed 
by Navarro et al. (1997), while u(x) = l/(x 3//2 + x 3 ) is suggested by Moore et al. (1999). Both the 
scale radius R s and the normalization 5 can be expressed in terms of a concentration parameter 
c(M) = R2oo/R s , where R200 is the virial radius within which the average halo density is 200 times 
the mean density of the universe. Specifically, R s oc M 1 / 3 /c, and 5 oc c 3 g(c), where g(c) is a 
logarithmic factor with g = l/[ln(l + c) — c/(l + c)] for Navarro et al. (1997) and g = 1/ ln(l + c 3 / 2 ) 
for Moore et al. (1999). High resolution simulations of individual halos show a roughly power-law 
c(M) with some, perhaps substantial, scatter (e.g., Cole & Lacey 1996; Tormen et al 1997; Navarro 
et al. 1997; Jing & Suto 2000), so we will write 

where M* is defined by a{M*) = 1, and (3$ is used to parameterize the uncertainty in the mass 
dependence, which has been found to be in the range < [3$ < 1/2. 

We now derive analytically the small-r behavior of equation (1), in which the factors dn/dM, 
R s , and 5 are all functions of the halo mass M. For small r or high k, the integral for £ih(r) or 
Pih(k) converges before the exponential cutoff in the mass function, and the dominant contribution 
to the integral comes from masses near the scale for which r/R s = 1 or kR s = 1. We then have 
P(k) w Pih(k) oc / dMv a u 2 {kR s )g 2 {c) at high k. Changing variables to y = kR s cx A;M /3o+1 / 3 , 
we obtain the asymptotic behavior 

ar)<xr~\ A(k)^Ank 3 P(k)^k\ ^ = ^ 2 (3^ + 1) ^ ' ^ 
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where we have {3 = (3o if the logarithmic factor <?(c) in 5 is ignored, or more accurately, we have 
(3 ~ 0.8 /3o if the effect of g(c) is approximated by a power law. From equation (4), we see that in 
order to reproduce the two-point stable clustering result for a given spectral index n, the indices a 
and (3 must satisfy 

18/3 - a (n + 3) _ 9 + 3n 

2(3/3 + 1) 5 + n ' ^' 

For n = —2, for example, this is satisfied if (3 = 0.25 or (3q ss 0.31 for a = 1, and /3 = 0.2 or 
A) ~ 0.25 for a = 0.4. 

For the three-point function £3 and the bispectrum B, the integrals in equation (1) converge 
if e < 3/2 for an inner halo density of r~ £ . For £3, B, and the commonly defined three-point 
amplitude Q oc B/P 2 , we find 

e 3 (r)ocr-™, fc 6 5(£;) oc A™ , Q(k) oc A; 73-27 , 73 = ^^±^ . (6) 

Stable clustering requires a scale-independent and 

a (n + 3) „ 

which holds only if a = or n = —3. Combining equations (5) and (7), we find that stable 
clustering requires a = and (3 = (3 + n)/6. All work on halo mass functions performed thus far 
has reported a much larger value for a; stable clustering therefore does not appear to be followed, 
which is consistent with the scale-dependent Q(k) reported in Ma & Fry (2000a). 

We illustrate these results by showing in Figure 1 the power spectrum and bispectrum com- 
puted from the analytic halo model for three sets of (a, (3) for the n = — 2 scale-free model. The 
values of a and (3 are all chosen to give the same slope 7 = 1 for the nonlinear two-point func- 
tion; the corresponding power spectra are therefore nearly identical at high k. The three-point 
Q, however, have very different high-Ar behavior for the three cases, and only the case a = pro- 
duces an approximately scale-independent Q. (Q is not perfectly flat for (a, (3) = (0, 1/6) because 
as mentioned above, the simple expressions in eqs. (5) and (7) are derived by approximating the 
logarithmic factor g(c) in the halo profile as a power law.) 

For the four-point and higher correlation functions, provided the integral is convergent without 
the exponential cutoff in the mass function, the same change of variable gives 

freer-™, **»->J*«*™. T,= 18 %^;« (n + 3> . (8) 
Stable clustering requires 

7K -( N -1) 7 = (J V-2)|L|1| = 0, (9) 

which again is satisfied by a = or n = —3. So if the three-point correlation follows equation (7), 
stable clustering continues to hold for all orders for which the integral converges. At large k, we are 
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probing the inner cores of halos, where the inner profile goes as r~ e . The Fourier amplitude u thus 
falls off as u(k) oc (kR s )~( 3 ~ e \ and convergence requires N[l — (3 — e)(3/3+l)/3] — l + a(n+3)/6 < 0. 
For e = 1 (Navarro et al.) and n > —2, this holds for all N, but for e = | (Moore et al.) and 
n = —2, it holds only as far as N = 3. 



3. Pairwise Peculiar Velocities 

Besides the scaling properties of the iV-point correlation functions discussed above, stable 
clustering also predicts that the mean peculiar velocity v between pairs of objects of separation 
r must cancel the Hubble expansion Hr on small scales so that bound halos maintain a fixed 
physical size. The velocity ratio v/Hr is related to the two-point correlation function £ by the pair 
conservation equation (Peebles 1980) 

-JL = j ^1 no) 

Hr 3(1+0 d lna ' y J 

where a is the expansion factor and £(x) = 3x~ 3 dx' x' 2 £(x') is the mean two-point correlation 
function interior to x. Stable clustering requires —v/Hr = 1. 

The analytic halo model of Ma &: Fry (2000b) makes specific predictions for the behavior of 
v/Hr. As eq. (10) indicates, the main task is to compute the time dependence of £. Equation (1) 
shows that the time dependence of £ enters through the rms density fluctuation a in dn/dM of 
eq. (2), and the halo concentration parameter c(M) of eq. (3). We obtain 

= / dM [f(u 2 - a) + (25' - 3 + A') c' ] ^ R* P X(r/R s ) , (11) 

where / = dlna/dlna s=s il a6 is the standard growth rate of the density field, 5' = dln5/dln c, 
A' = dlnA/dln x, and d = dlnc/<ilna. Specifically, for the Navarro et al. profile, 5(c) = 
(200c 3 /3)/[ln(l + c]-c/(l + c)] and X(x) = 8ir[(x 2 + 2x + 2) ln(l + x)/(x 2 + 2x) - l}/(x 3 + 2x 2 ); and 
similar expressions for the Moore et al. profile can be found in Ma & Fry (2000b). The term d for 
the time dependence of the concentration parameter is somewhat uncertain. Numerical simulation 
results of Bullock et al. (1999) give a stronger redshift evolution, c cx a, than Navarro et al. (1997), 
while more complicated forms have also been proposed (e.g., Cooray et al. 2000). We note that in 
equation (3), c(M) has an implicit dependence on a through M* which roughly agrees with Bullock 
et al. when (3 = (3 + n)/6. 

The prediction of the halo model for v/Hr can be calculated from equations (1) and (11). 
Before doing so, it is useful to derive first the asymptotic behavior of v/Hr in the deeply nonlinear 
regime at small r. On such small scales, £ 3> 1, and equation (4) gives a power law £ oc r 1 and £ = 
3£/(3 — 7). For a scale-free model, the self similarity solution gives d^/dlna = — 2/(3+n) d^/dlnr. 
We find that equation (10) then reaches 

v rwo 2 18/3-a(n + 3) 



Hr ra + 3 6 + a (n + 3) 



(12) 
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It can be easily checked that —v/Hr approaches unity, as it should, for any pair (a, (3) that satisfies 
the stable clustering condition of equation (5) for the two-point function. 

Figure 2 shows our results for the general behavior of —v/Hr calculated from the analytic halo 
model. An n = —2 scale-free power spectrum is used for illustration. Although £ is dominated by 
the 1-halo term £ih of equation (1) at £ > 1, which is the range of interest for this paper, we improve 
the accuracy at £ < 1 by adding the linear theory £j oc r ~( 3+n ) oc r _1 to £1/4. As Figure 2 shows, 
the resulting v/Hr indeed follows the linear perturbative prediction at £ < 1, and equation (12) 
at £ 3> 1. We find —v/Hr — > 1 for the same pairs of (a, (3) shown in Figure 1, which were chosen 
to illustrate the stable clustering scaling of the two-point function. However, —v/Hr reaches a 
smaller value at £ 3> 1 for the (a, (3) suggested in the literature, e.g., a = 1 or 0.4 for the mass 
function (Press k Schechter 1974; Sheth k Tormen 1997; Jenkins et al. 2000) and (3 = 1/6 (Cole 
k Lacey 1996; Navarro et al. 1997), implying that the pairwise peculiar velocities are too small 
to cancel exactly the Hubble expansion. As a consistency check, we have computed the derivative 
d^/da using both eq. (11) and [£(ct + Aa) — £(a)]/Aa. The two methods give identical results, but 
eq. (11) offers a more physical insight into the various contributing terms. For example, the broad 
peak in —v/Hr for 1 < £ < 100 is the result of the term in u 2 in eq. (11), which reflects that the 
integral on these scales is dominated by mass scales near M* where a(M*) = 1. 



4. Summary and Discussion 

In this Letter, we have derived the asymptotic nonlinear behavior of the iV-point correlation 
functions £/v (eqs. [4], [6], [8]) and the pairwise peculiar velocities —v/Hr (eq. [12]) in the framework 
of the recently proposed analytic halo model of Ma k Fry (2000b). We have shown that their small 
scale behavior is consistent with the stable clustering hypothesis only if dark matter halos satisfy 
certain criteria. The two halo parameters whose variation we have explored are the logarithmic 
slope a of the halo mass distribution dn/dM at the low mass end in eq. (2) and the slope (3q for 
the mass dependence of the halo concentration parameter c(M) in eq. (3) (recall (3 ~ 0.8(3o). These 
two parameters are highlighted because results to date from numerical simulations have indicated 
a significant uncertainty, with 0.4 < a < 1 and < (3 < 1/2 being plausible values. 

From the derived asymptotic nonlinear behavior, we have obtained analytically the relations 
that must be satisfied by a and (3 in order for stable clustering to occur. Equations (5), (7), and 
(9) summarize the results for the iV-point correlation function £at. These equations and Figure 1 
indicate that although the two-point stable clustering condition alone is satisfied by an infinite set 
of (a, (3) for a given n, the three- and higher-point stable clustering conditions, £jv oc ^ N ^ 1 , N > 
3, imply a must be zero. Achieving stable clustering to all orders therefore requires stringent 
conditions: a = and (3 = (3 + n)/6. For non-scale-free models (e.g., cold dark matter), n is the 
index at the high-fc end of the power spectrum. 

Compared with the hierarchical model of Davis k Peebles (1977), results from this paper and 
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Ma & Pry (2000b) show that the analytic halo model makes more general and physical predictions 
for the behavior of the iV-point correlation functions and the pairwise velocities. Imposing the 
stability condition a = and (3 = (3 + ra)/6 in equations (4), (6), and (8), we indeed recover their 
results of £ oc r 7 with 7 = (9 + 3n)/(5 + n) at small r and £at oc £ Ar_1 . The condition f3 = (3 + n)/6 
can be understood further by examining individual halos. Using c oc Mf°, M* oc a 6 /( n+3 ), and 
/?o ~ /? up to a logarithmic factor, we see that (3 = (3 + n)/6 gives approximately c oc a(t), which 
is consistent with the simulation results in Bullock et al. (2000). The halo scale radius R s is then 
approximately constant in time in physical coordinates, as expected for a stable system, and the 
physical density p of individual halos is indeed very close to constant in time, more so for very 
small halos that form early (;/ > 1, c > 1) than for those that form recently {y ~ 1, c ~ 1). 

Yano & Gouda (1999) have found that for halo density p oc r~ e with e < 3/2 (such as Navarro 
et al. and Moore et al.), the velocity parameter v/Hr approaches and the two-point correlation 
function £ approaches a constant (i.e. 7 = 0) in the nonlinear limit. Neither behavior is seen in 
our Figures because their derivation ignores the mass distribution function dn/dM and is therefore 
valid only for equal mass halos. 

Our earlier paper (Ma & Fry 2000a) has already shown signs of departure from the stable 
clustering hypothesis in high resolution iV-body simulations. In this paper, the specific values 
a = and (3 = (3 + n)/6 required for stable clustering provide additional evidence that this long- 
cherished hypothesis may not be applicable in all situations. In such cases, one consequence is that 
the frequently used fitting formulas for the nonlinear power spectrum (Hamilton et al 1991; Jain et 
al. 1995; Peacock & Dodds 1996; Ma 1998) will need modifications at high k, e.g., k/k n i > 50 for 
the n = — 2 model and k > 20 h Mpc -1 for cold dark matter models with a cosmological constant 
(see Figs. 3 and 4 of Ma & Fry 2000b). We believe the analytic halo model is a powerful tool that is 
providing new insight into the nonlinear regime of gravitational clustering, the most fundamental 
process in cosmology. 

C.-P. M thanks Jim Peebles and Uros Seljak for useful discussions. She acknowledges sup- 
port of an Alfred P. Sloan Foundation Fellowship, a Cottrell Scholars Award from the Research 
Corporation, a Penn Research Foundation Award, and NSF grant AST 9973461. 
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Fig. 1. — The power spectrum (upper panel) and three-point amplitude (lower panel) given by the 
analytic halo model of Ma & Fry (2000b) for the n = — 2 scale-free model. Three sets of (a, (3) are 
shown, where a and (3 parameterize the halo mass function and halo concentration as defined in 
eqs. (2) and (3). (Note (5q in eq. [3] is related to (3 here by /3 ~ 0.8/?o; see text.) The Navarro et 
al. profile with c(M*) = 9 is assumed. Dotted straight lines show the lowest order perturbative 
predictions. While different combinations of (a, (3) give nearly identical A(fc) which all agree with 
the two-point stable clustering result (see eq. [5]), the three-point Q(k) has distinct high-A; behavior 
(see eq. [6]), and stable clustering holds to all orders only if a = and (3 = (3 + n)/6. 
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Fig. 2. — Prediction of the analytic halo model of Ma & Fry (2000b) for the ratio of pairwise 
peculiar velocity and Hubble flow, —v/Hr, vs. the two-point correlation function £ for the n = — 2 
model. The velocity ratio approaches a constant in the deeply nonlinear regime (see eq. [12]), but 
the value is not always unity as required by stable clustering. For example, —v/Hr — > 1 at £ 3> 1 
for (a,P) = (1, 1/4), (0.4, 1/5), (0, 1/6), but it reaches a smaller value for a = 0.4 or 1 and (3 = 1/6, 
which are the range of values found in simulations. Dotted straight line shows the linear theory 
prediction, which is followed accurately at £ < 1. 



